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numerical Method for Solution of Systems of non-Stationary Spatially 
One -Dimensional Nonlinear Differential Equations 


S. K. Morozov and 0. P. Krasitskiy 


We propose a computational scheme and a standard program for solving 
systems of nonstationary spati^ly one-dimensional nonlinear differential 
equations using Newton’s method. The proposed scheme is valuable because 
of its universality and the fact that it reduces to a minimum the work of 
programming. We give a detailed description and present the text of the 
standard program which realizes this computational scheme. The program 
is written in the FORTRAN language and can be used ■without change on 
electronic conputers of type YeS and BE3^6. The standard program 
described permits us to find nonstationary (or stationary) solutions 
to systems of spatially one-dimensional nonlinear (or linear) partial 
differential equations* The proposed method may be used to solve a 
series of geophysical problems which take chemical reactions^ diffusion, 
and heat conductivity into acco^mtJ to evaluate nonstationary thermal 
fields in twu-dimensional structures when in one of the geometrical 
directions it can take a small number of discrete levels, and to solve 
problems in nonstationary gas d 3 mamics. 


INTRODUCTION 

In contemporary science and engineering practice, one often encounters 
situations when a specialist not familiar with the methods of computational 
mathematics and the fine points of programning must solve on an electronic 
conputer one or another mathematical problem. In this case, he turns to 
a libraiy of standard programs. For example, standard programs have 
yielded a great extension of the solution to the Cauchy problem for systems 
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of ordinary differential equations by methods of Runge-Kutta type* The 
level of the developments of present-day computational engineering permits 
us to pose a question about the creation of standard programs for more 
complex mathematical problems* In this work, we propose a standard program 
for solving the boundary problem for a system of nonet at ionary spatially 
one-dimensional equations* A very wide class of physical problems can 
be reduced to such a ^stem* For exan^jle, many models of geopl^sical 
phenomena, in which diffusion (heat conductivity) is the dominant factor, 
are dften formulated mathematically as nixed problems for equations of 
parabolic types either initial or boundary conditions are given* As a rule, 
the ^Inknowns are related to each other in a nonlinear manner, which signi- 
ficantly complicates the solution* 

At the basis of the method is the use of the implicit difference scheme* 
To solve nonlinear difference equations at each time step, we use Newton* s 
iterative method* At each iteration, a linear system of algebraic equations 
is solved by a modified method of Gauss, taking into account the sparseness 
of the.rimatrix, by choosing a pivot element by column and a normalization 
by row* They ^behave :in ^ analogous wsy in problems of mathematical 
chemistry (see collections Ul ,2 ] !) » In gephysics a similar method was used 
in work |]3jj* 

The ^stem of equations under consideration can be very stiff, :i-*e* 
it can contain time constants with essentially different values. For 
example, in geophysical and astrophysical problems describing the distribution 
of atmospheric gaseous components by height (the determining processes are 
chemical reactions and vertical diffusion), the stiffness depends on the 
great difference in the rates of chemical decay of the various components, 
and also on the coefficients of diffusion at various heights. It is known 
ill] .[that, on account of a stiff restriction intrinsic to them at a step of 
time, it is unreasonable to use explicit difference schemes to solve stiff 
systems of differential equations, because the utilization of such schemes 
requires a large outlay of machine time while, in most cases, such ^sterns 
in practice cannot be solved by explicit methods on present-day electronic 
computers* 

The proposed scheme is especially advantageous when the determination 
of an exact solution in time is not required, but we are required to 
find the steady-state behavior* In cases of the lack the * smallness ■ F] 
of the coefficients of the spatial derivatives, as this often occurs in 


2 



ORIGINAL PAGE IS 
OF POOR QUALiry 


problems of mathematical chemistry, the scheme can be shown to be unstable 
(the scheme is A-stable * In such cases, we need, to watch for negativity 
of the eigenvalues or, at any rate, for positivity of the solution, and to 
control the pace of time or to use further iterations, excluding the result 
of the solution in the negative domain (see, for example, the article of 
Yu»M* Volin et al. and the article of B.V. Pavlov in collection! [i].) 


Linearization by Newton* s method permits us to provide an effective 
iteration process for solving nonlinear problems, to easily guarantee the 
conservativity of the difference scheme (a property of the scheme without 
which it is often impossible to obtain a numerical solution to the original 
system of differential equations). 

It is well-known how much work and time it takes to write down and 
debug a program. Therefore, by writing dom the program, we set up the 
problem, reducing to a minimum the work of programming and ensuring the 
possibility of an operational introduction of changes in the original model 
(for example, in geophysical problems the addition of new components and 
chemical reactions). This problem is solved by the calculation of a Jacobian 
matrix by means of numerical differentiation. In this way, we decrease the 
probability of error in the written program, while the user is freed from 
tiresome work and can concentrate, his attention on the physical formulation 
of the problem and on the numbers of the e^qjeriment. 




1. THE MATHEMATICAL METHOD 


Let there be given the system of equations 


( 1 ) 


where k,i = 1, ...,M (i.e. into f^^ enter terms ^idth different i = 1,...,M), 
M is the number of equations, yj^ is an unknown function, t is time, and x 
is the spatial coordinate. 

We are given some initial conditions 


2) The definition of A^stablity can be found, for exanple, in the article 
by L.S. Polak et al. in the collection [2 ]> 
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( 2 ) 


(0,x) ^ Pk Cx) , . 

We seek a solution on the interval a < x < h* /6 

•M mm mm 

The boundary conditions have the form! 


^ (a) , 


(a) X 

dx J^o 

dx J — 0 


> 


> 


( 3 ) 

( 4 ) 


V7here i,k = 1,...,M. 

We pass from continuous variables to discrete ones. For this we choose 
a uniform mesh in x* 


- aHn-D-tX ; ' c^X^d>-a)/W-0 ;■ ( 5 ) 


(if it is necessary to introduce a nonuniform mesh, xre can usually make a 

substitution of x by x*, so the mesh in x’ will be uniform). The time 

steps will be numbered by means of the index j, vbile the moments of time 

coi*req)onding to the steps will be denoted t.. By means of y^ we denote 

3 K,n 


approximate values of the continuous variable y, (t,x) at the point (t.,x ). 

The derivatives entering into f^^ are approximated by difference relations 
of not more than three points x^. For examples 



( 6 ) 

( 7 ) 


while the derivatives entering into and ^ are approximated by the 

relations 


Ik- (hi, a) 'jL-yji. 

dx / • ~ ^ 3 : 


( 8 ) 


k 
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/rr 


X 


( 9 ) 


In order to avoid restrictions on the time step ^ ~l> connected td.th 
instability^ vre make use of a two-tiered implicit schemes 

yLn^ ’'^in - fuj . Ui 

^ V • — "" ■ ■ 

*ere is an approximate value of f^^ at the point n at moment of time 
to which is obtained by the substitution of expressions(6) and (7) into 
f for n = 2, . . . ,H-1. 

The botmdary conditions can be approximated by the expressions 

"^,Uyi^,yiz)=o , (n) 

, c«) 

where and are approximate expressions for and ' ^ , 

which are used in the substitution of espressions (S)j(9) into (3)i(4)» 
Equation (lO) has first-order accuracy in time* To increase 

the accura<^ in time# we can put the right-hand side in the form 




(13) 


forms 


d^i jt;,a) J -y j. 

3x 7 Z AX ^ 


(14) 


t? 


For p = 0*5a equation (10) ^pixDxlmates the (Jifferential equation (l) with 
double precision in time. For g <0*5, the difference scheme becomes 
unstable. 

To increase the accuracy of approximating the boundary conditionsj 
the derivatives entering into ' and can be written in the 
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2a X 


Cl5) 


In this case the deriiratiTes mentioned above are approximated with double 
precision in x. 

Eq,uations ClO) , Cu) » and {.12} form a nonlinear system of algebraic 
equations 5 consisting of N*-13 equations and N*M unknowns. We will solve this 
system by lewton's method, i.e., -at each, iteration S we solve the following 
system of linear algebraic equations 

'<>. 5 ^/ uH 


/8 


U ' u< 


4 r; 


M h 






y y 

l_ Z- 3u/,5 

,-_y O^n . . . - 


Ui . f=/. - 




(16) 



; w M 




(17) 


A/ 








K 


' ££ 3y{^, 








(18) 


For S = ‘0, at the zero-th' approximation the values begin with the previous 

U 


time period: 


The derivatives of the functions f^^ ,n- , and 


with respect to the unknowns (of the Jacobi matrix) are found by numeri- 
cal differentiation. 

The -, nume r i c al differentiation is effected by the formula 

For clarity and convenience of our further presentation , we represent 
the system .of equations (l6) , .(17), and (l8) in matrix form: 


\ fifj &f/ ^ ‘ ‘ ^ ^ 0 i 0 


F, 


Oa/ 

j S//-/ CtH 0 • ’ . • 0 0 0 I 0 


F,.i 



' 0 bfj.% • 0 0 0 0 

r * * * 

• *^ * • * • * ‘ ^ , 

* • * 

• 

Fn-z 

• 



• ’ • • . ♦ • • • 

1 

J 

. • ♦ • • « # 

w • ^ • • • 

0 0 0 0'* . • 0 


F^- 


• 

% 

0 0 0 0’ 





0 0 0 0.* 0 Oj 


F^ 




( 19 ') 
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Here A, G, and C are squaxe matrices of size and F and D are M-dimensional 

vectors. In this form, the matrix has a "block.— tri diagonal form under the condi- 
tion that each- of its "block elements, in turn, is a square matris. In the upper 
and lower rows 'Of the matrix eq.uation (l9) we write eq.uations (l8) and (l7), 
while in the interior rows is eq.uation Cl6). The vector 

Expressions for the elements of the matrices A, B, and C and the compo- 
nents of the vector D have the form: 


1 ) for n = IT and i ,k = 1 , . . . ,!M 


3 






hf tTi 


■ .. . i,A'+£-3 

2 } for n = 2 , . . . ,lT-2 and i ,k = 1 , . . . ,M 

i - yi'l. V f" 


C -1^ 

> ^ yhi 

UA'. 


K,n. r — U <^r:.rv 

us ’ ayl^ 

,,rh 


= 


i-./ C-f . 


where 6. . is the Krone cher delta symbol; 

3) for n = 1 and i ,k = 1 , . . . ,M 

-'n ' .^zZjl X 


c -- ^ 

3y('i ■ > 


d t V"'" ■ 

k, t"/ 

The linear system i(l9) at each iteration S is solved by a modification 
of the. method Gauss, which takes into account the tridiagonality of the matrix, 
with the choice of a pivot element by row and a normalization by column [5] • The 
basic idea of the elimination method consists of the fact that we use a linear 
transformation to make the diagonal elements of the 
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matrix into ones, and the below- diagonal elements zeros. This 
is done sequentially from top to bottom. To demonstrate the 
recurrence formula for passing from n to n-1, we introduce 
the rectangulsir matrix 





The meaning of the matrices B, C and the vector D will become 
clear from further presentations. 


We take the first column of the matrix B2, find in it the 
maximal element as a model, and place the row in which' it is 
found in the first position of the matrix, BZ. We divide each 
element of this row by its first element. Now, subtracting the 
first from the other rows , after multiplying by an appropriate 
factor, we obtain zeros in the first column of these rows. We 
find the maximal element in the second colxmn among all rows 
except the first , and by means of this row xre obtain zeros in 
the second coltimn of all rows except the first. We conduct an 
analogous transformation, further, vTith the third column, etc. 
As a result the matrix BZ is reduced to the form 



E 

■a 

Qh 



0 


Cn^ 




1 


where E is the identity matrix. 
For n = N, \'7e have 



1 

S;/ 

Cf/ 

1 


1 



j 


and the passage from N to N-1 is included in the standard scheme 111 
of the passage from n to n-1. 

For n = 1, before the standard scheme of passage is applied, 
is eliminated by means of the identity matrix obtained for 
n = 2, In other words, the matrix 
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by means of the matrix 




||£ P. ^ 


reduces to the form 


-11 


|o e>: 


Then the standard scheme is applied to the matrix 


B 2 = 


Bs^ C2 0 2)5 


b; c;o^ 


As a result, we obtain 


- 


E 

\o t; g s; 


from which 


and then 




p^yR^-pz-F., . 

And finally the recurrence formula for finding the solution has the 


form 


Qa'f^n-2 


Thus we find an approximate solution to the boundary value problem 
(1-4) at the moment of time t. in S+i iterations. If the iterations 
correspond to a given precision, then we can continue iterating equation 
( 1 ) in time. 


2. DESCRIPTION OF THE PROGSiAM 


The program permits us to find a nonstationary solution y^(t,x) and 
a steady-state, solution y^(x). The nonstationary solution can be found by 
integrating equation (l) with a constant step in t or with an automatic 
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choice of step. !Di the case of an "automatic ehoice of step"^ the time 
step is chosen depending on the number of iterations with reqpect to 
nonlinearity at the given step. Thusj the controls on precision do not 
depend on time* 

The steady-state solution can be found by the method of settling down, 
i*e. we integrate in time up to. the moment of time exceeding the maximal 
characteristic time of the problem (better with an automatic choice of 
step)* or by iterations without settling down in time. In the latter 
case* we can solve directly a system of stationary equations* i.e. 
equations (l) without time derivatives. 

Intobthe standard part of the program* we introduce three subroutines! 
HOK* OUT* andlSCL. 


The sub37outine HOK controls the iterations* the choice of .step in 
time*'.VJ the printing of the result* etc* 

The subroutine OUT effects the printing of the result. 

I calculaten the coefficients for the imknovms 


In the subroutine ISCL 


we 


in "A linear system of difference equations* and this system is solved by 
Gauss" method. 

The user writes the subroutine FF* where the difference equations are 
witteni^.down* and the basic program* where values are given to the physical 
constants and from which originates the access to the standard part of the 
program. 

The' access to the subroutine has the form* . 


Here Kf is the number of partition points of the x-interval. 

M is the number of equations. 

IB. is the number of steps TAU in time* across which the printout of 
of the solution is produced. When it is necessary to print out 
the results at fixed moments of time t^* xre need to sequentially 
access the subroutine HOK* giving TKOH = 

ITAU — 0 is the signal of a calculation with constant step TAU in time. 
ITAU 7 ^ 0 is the signal for a calculation with automatic choice of time 

step. The step is increased by 1.5 times when the number [of iter-; j 
ations at the previous step equals 1 and decreases by 2 times 
when it is greater than 3* Thus* there is no linear relationship 
between the choice of step and the accuracy of the solution in time 
When we need to obtain a steady-state solution by the method of 
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settling down, then we should set ITATJ ^ 0 . 

ITER is the maximum allowable number of iterations at the step TAU^ 
when the number of iterations IT exceeds ITER, there occurs 
a decrease in the step TAU by two times. 

ITB= 1 is the signal to print the solution at each iteration. 

EPS is the relative accuracy of the convergence of the iterations. The 
magnitude of the increment^"” 'in the numerical differentiation 


is given by- 'the formula 


- 0 (for example, when the initial approximation!. 


When Ay. 

i,n 


t0- = O) y. = EPS ( in a program with double precision Ay. \ 

3 . 3 . ^ 13 . 

= (EPS)2). _ ^ 

TKOH is the end of the interval of integration in time; when TAU > 

TKOH/iO, the step TAU does not increase; when TAU < TKOH/iO^, 
the calculation ceases and the statement "THE TAU IS 3IALLER 
THM ALLOWED" is printed, and the values of the current moment 
of time and TAU are given. 

TAU is the value of the time step! the initial step or’ constant step, 

depending on the value of the parameter ITAU. V/hen TAU = 0 , yie 
seek a steady-state solution of the system directly by iterations 
without a settling-down in time. In this case, when the number 
of iterations IT > ITER, the statement "IT = ITER ITERATIONS DO 
NOT CONVMGE VIETH THE GIVM ACCURACY EPS" is printed and the 
values of the next two iterations are given. 

Y0 is the dimensionality file W-M of the initial conditions ;-i 

fbr„ the output of the subroutine HOK,Y0 is a solution for 
t = TKOH or in case TAU = 0 .i* a steady-state solution, 
is the left enc^oint of the interval of x^ s. 

XL is the right endpoint of the interval of x* s. 

The subroutine HOK works with the subroutines ISCL, OUT, FF, 
and the library subroutine ABS. 

ISCL is the subroutine for solving the system of linear algebraic 
equations ,. An accidental result of the form "MATRIX IS 

DEGENERATE....... " is possible. Such a result, as a 

rule, indicates an ingsroper formulation of the problem (an 
'incorrect assignment of initial or boundary conditions) or 
errors in the subroutine FF (improperly written equations, errors 
in the coefficients, etc.) 
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OUT is the subroutine of printing the solution. The value of the current 
moment of time is printed BP = ...J the column of x v^ues and 
the columns of values of the unknoim functions yj^ are also printed. 

FF is the subroutine which is written by the user. In it are given the 
formulas for computing the right-hand side of f^ and also of 

and • In this subroutine we must describe the 

general block HIFs 

CQMMOh//HlF/BP,H, TfiVi, F ( iO), 

- - - 


where BP is the moment of time t^| H is the coordinate step^ 
TAUl is the reciprocal value of the time step TAU^ XH‘, 


.is the file of x • for n-= 1,. 
. ‘ n 


Y0 is the value of 


the function at level t . , } and Y is the value of the function 

J— X 

at level t^ of the previous iteration. 


Access to the subroutine FF originates from the standard 
part by means of the operator GALL FF(K>NN)j where W is the 
current index ' of the point x^, M is the collection of partition 
points N Correspondingly^ the first operator of 

the subroutine FF must have the forms SUBROUTINE FF(N,M)> where 
M, NN are formal parameters having the sboveMnentioned meanings. 

The connection between the subroutine FF and the basic program 
can be effected through additional general blocks (also inclucJing 
• unmarked ones } . ' • ; i \ 

On the basis of the data, we calculate the 
of F of dimension M, the elements of v/hich are given by the 'flle.| 
numerical values of expressions representing themselves or , ellther.' 
difference approximations of the right-hand sides of the 
equations at the point x^, or a difference approximation of the 
boundary conditions. 

The variables entering into the general block HIF cannot 
be used in the subroutine FF as identifiers of other variables. 

To economize machine time, it follows that we should avoid 
in the subroutine FF calculations of expressions not depending 
on j and n, because to calculate the Jacobian matidx at each 
iteration at each point x^ is the result of M+1 accesses to the 
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subroutine FF,< ..Purtheiwrej if these expressions enter the sub- 
routine PP, then they also are computed M4-1 times. 

In the operators of description COMMON and DIMENSION, the | fiies.< 
of the variables are written so that we can use the program when N < 21 and 
M < 10. When it is necessary to increase N or M, it follows that we should 
insert corresponding changes into these operators. In different programs 
the operators have the following forms 

COMMOM/Hir/ 6P,H, T/)U XPffl/) , Fffi ) , 

in the subroutines HOK, ISCL, PF and OUT; 


2) DIMENSION YJ^(N,M) in the basic program; 

3) DIMENSION Y^(N,M) in the subroutine HOK; 

4) hlMEMSIOA' 52 ?2f^, P2 ('/V-/,/Y, /i), ■' 

, 

'in the subroutine ISCL. 


Here the variables N and M and expressions containing them need to be replaced 
by appropriate numbers because the sizes of the Idimensions of the fil e's are ',j 
given by integer constants without sign. 

For example, if N == 51 and M = 5> the the operator DIMENSION in the 
subroutine ISCL needs to be replaced by the operator 


mMeNSIOi^ 3ZfrO/6l Pl(50,S,S\ QZ (vs,5,5h CZ&, Pi I 5 ) . ; 

We need to make analogous changes in operators 1-3* After we have made these 
changes, we can use the program for arbitrary N < 51 and M ^ 5* 

The necessity of changing theropera'^rs'.'of description COMMON and 
DIMENSION for increases^in N or M creates definite inconveniences in the 
work and is related to the limitations of the programming language Fortran, 
in which, in contrast to programming in Algol, we cannot use dynamic 
/files . To eliminate this deficiency and also to fi^rease P 
the econony of the programat.,-we should transfer from multidimensional 
'files , to one-dimensionaL ones”, for example is done in the ^ ^ 

IBM library of standard programs {6}. The authors hope to return to this 
question at some time in the immediate future. 
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Example 


Let there be given a system of linear differential equations of parar- 
bolic type* 


3 SO 

_ 9 V 

3 t 

”” 2. 
O vC 

39 z 

_ 9 % 

3t 

dx^ 


ri- y, 


2 .? 


( 20 ) 


where 0»5 ^ x < Ij t > 0. 

We put on this system of equations the following boundary and initial 
conditions t 


d.y^f (t, 0 ,d) Jq 


O’X 


> 


5 X / “ ^ 

V/ io,%) = 0 


y^(t, i)-o , 


( 21 ) 


The analytic"': solution to the problem formulated has the form* 

9f(t,x)= t (%x) , 


7 o 

^ 6in (nx) . 


( 22 ) 


We choose a uniform grid in x and ts 

; c^x-0,0Z5 ; 

^00- i^t=0,00-f . 

The differential equations (20) may be approximated by the difference 
relations? 
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The difference expressions for the boundary and initial conditions 
(21) have the form 


x-3 yfj V yjs. ~ 

2 A X 

2ax , 


=o v 


^'0 > 


; i4%=t3mr5i-Xj 


where j = 0, « . . , 100, n — 1, . . . , 21. 

l^en we accessed the standard program HOK, the following actual parameters 
were - given s 

ti = 2x7 7 tV^ 18= 50, ' if 

EPS- o.ooi. TKOH=.o,i,_ TW= 0 , 001 , , o,5, W • i,; 

In the subroutine FF, the right-hand sides of equations (23) and the 

bo^lndary conditions (24) are widtten down* In the appendix we present the 

the text of the subroutine FF, where H = Ax, H2 = (Ax)®, Y(k',N) = y£ , 

*1 X ' n 

;Y0(K,n) = N = n. N = 1 corresponds to x = x0, and N = UN corresponds 

to X - xl, 1 < N < NN. 

The initial conditions are assigned in the basic program (see the text 
of the basic program in the appendix). 

In the table we cite the epproximate and exact solutions at the moment 
of time t = 0.1. 
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Solution 


Coordinate 


0. 50000 

0,52500 

0,55000 

0,57500 

0 , 6C000 

0,62500 

0,65000 

0,67500 

0,70000^ 

0,72500 

0, 75000 


0,87500 


0,95000 


1,00000 


Approximate 



Solution jTg 


0,03726^ 0,037271 

0,037149 0,037156 

•0,036805 0,036812 

0.036235- . 0,036241 

0,035441 . 0,035447 

•0,034428 • 0,034434 

0,033203 0,033209 

0,031773- 0,031779 


0,030148 

0,028396 


0,016918' 

0,014261 


0,030153 

0,028341 


0,026350 0,026355 


0,77500 . 0,024202 0,024206 

0,80000 0,021904 0,02X907.. 


0.0I947T 0,019474 


0,016921, 

0,014263* 


0,011515 ■ 0,011517 
0,008699 ; 0,008701 
0,005629 0,005831 

0,002924 0,002924 


0, 37300 
0, 37186 
0,36842 
0, 36271 


0 . 3^^62 

0,.:^3236 

■0, :|I605 

0, jpi78 
< 

f), 28365 
0,26376 
24226 

L 

vu ^ ^ 



0,02924 











Before solving conplicated problems by the standard program described 
here, we recommend to the user that he solve one or more simple problem s, 
the solutions are which are well-known* 



COKCLUSION 

The method described here was shown to be effective in solving stiff 

systems of differential equations for modeling the chemical conposition 

of the atmosphere of Mars» i [Tl^ l^stem|n)]of equations of parabolic type 

was'- solved* An. iitplicit conservative difference scheme of second-order 

accuracy in \ Ax'-h ^as used* Steady-state and nonstationary solutions were 

found. The characteristic times of chemical decay for different components 

in this problem varied between the limits 10“ to 10 see*, while the 

characteristic times for the establishment of diffusion equilibrium varied 

2 7 

depending on . ‘ height within the limits 10 to 10 sec* In calculations 
with usual accuracy on electronic computers of type YeS and BSSIt-6, we obtain 
larger roundoff errors for time steps /kt of 10^ to 10 sec. These errors 
can be e::q)lained by poor conditioning on the system of equations (19)^ which 
is a consequence of the great difference in characteristic times of the pro- 
blem. It turns out that a steady-state solution is obtained only when 
calculating to double- precision . In this case there are no restrictions 
on the time step because the corresponding linear problem (19), which 

is obtained as a result of linearizing the nonlinear difference equations, is 
stable. ¥e consider the steady-state solution to be lobtjained when t 10^^ 
sec* (the maximal characteristic time of the problem ~ 10 sec*). 

Using the given method to solve the described system of equations ^21 

permits' us to obtain a series of nevr results on the chemical con^josition 
of the atmosphere of Mars. 


Pbr an abrupt change in the solution at a given time At, the emergence 
of negati'VB concentrations and instablity is possible* In this case, the 
iterations do not converge, and the program automatically changes the time 
step xmtil the concentrations become positive. 
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The proposed method is adaptable to the solution of problems 
of nonstationary gas dynamics. For example, to solve two-dimen- 
sional problems of the atmospheric heat circulation of planets 
[]8,9,10], an implicit scheme was used in which the angular 
derivatives were taken with preceding iterations,- However, to 
solve the one -dimensional boundary problems obtained from it 
along vertical lines, a' variant of the i^ro posed method in which 
the Jacobian matrisg was calculated anal 3 rtically was used. This 
method is also adaptable to the calculation of thermal fields 
in two-dimensional structures and other problems. 

In conclusion, v?e enumerate the basic merits and defects of 
the proposed method, and v/e give recommendations about its use. 
Regarding the merits of the method, we can list 

1, The possiblity of using the i>rogram for solving -^a 'Wide 
variety of physical problems , among which are those 
having a large range of characteristic tumes, 

2, Rapid convergence of the iterations for nonlinearity, 
owing to the use of NeT-jton's method. 

3, The possibility of easily constructing a conservative 
difference scheme, 

41 The automatic linearization of nonlinear difference 
equations by means of numerical differentiation, vj-hich 
significantly reduces the work of programming. 

It follows that we should consider as an essential defect 
of the method the fact that the calculational time (the number 
of arithmetic operations) in solving a system of M equations 

3 

is proportional to M , It follov7S that we should stress that 
this defect is related to the use of Gauss* method, i.e. it 
turns out to be closely interrelated with the merits of our 
method , 

Owing to, on the one hand, the stability and rapidity of 
convergence of the iterations, and on the other hand ' the 
reduction to a minimum of the programming v/ork of the described 
method, it can find xTide application in scientific and engineering 
computations . We note that more economical methods , in which the 
calculational time is i^roportional to , for example, often 
possess slower rates of co&^rgence of the iterations and, as a 
result, can prove less economical. 
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Recoromendatlons for Us& 0 P pOOft 

1 . Linear equations , 
a« Stationary equations . 

To solve such problems, it is necessary to solve only once 
a corresponding systemiJof difference equations, i»e« we need one 
iteration. For that one iteration to be made, it follows that 
we should set ITER = 1. However, by convention, it is useful 
to set ITER > 1 (for example, ITER = 5) and ITB. = 1, Then when 
^there' are no.' errors.^' -one superfluous iteration mil be made. If, 
however, there is an error, then with a significant probability, 
the iterations will not converge. The result of each iteration 
will be presented in print . ' ' ' ‘ , 

b. Nonstationary equations. 

Because there is in the program no control over the time 
precision, then in 'order, to guarantee it, it follows that we 
should solve the problem several times with different constant 
time steps \ At [ and should equalize the results obtained. 




2. Nonlinear equations, 
a. .Stationary equations. 

If the equations are strongly nonlinear and the form of 
the solution is unknovai beforehand, then it is necessary to 


solve the problem by the method- of bisection. In other vrords , 
it is necessary to give an arbitrary initial approximation and 
time step [_A^I, V7hich must be smaller than the minimal charact- 
' erisitic time of the problem. We need to undertake the calculex.' 
ation with an automatic choice of step up to a moment of time 
larger, j^han the maximal characteristic time of the problem i" 


b. Nonstationary 



Here the same, remarks as those in item "b" for linear 


equations are valid. 


/23 


We can use the method to solve multidimensional problen^ 
by methods of coordinate -wise decomposition. 

In the case of solving systems of equations having a large 

4 ' 

difference in characteristic times for different fxmctions, i.e, 
stiff systems, we srecomraend the use of a program with double 
precision, the text of which can be found in the a^^pendix. -It 
follows , in view of this , that such a program will occui>y twice 
as much space in the machine's memory. On an electronic comi)Uter 
of type YeS, the calculation time for such a program cannot be 


i§ 



practically increased, but on a B‘ESM-6, the space is increased 
five-fold « On the BESM-6, a number v/ith ordinary precision has 
12 decimal places, but on a computer of type YeS only 7. There- 
fore, for calculations on a computer of type YeS, we always 
recommend the use of double precision, because, on account of 
the restricted number of spaces , round-off errors which are too 
large can occur in numerical differentiation and in solving 
a system of algebraic equations , 

To solve systems of differential equations ^^th derivatives 
no higher than first order as part of the boundary condJitions, 
we can use difference equations written on one of the boundaries « 

In this case, the difference analogue of the time derivative 
is written do^m in a corresponding expression for or 

^ (see formulas (ll),(l2)). 

If an equation not containing a time derivative enters the 
system, then it folloTvs that we should multiply this equation 
by a sufficiently small number. 

We recall that the choice of a difference scheme belongs /24 
to the user. We have only the follovTing requirements for it: 

a) a scheme of ordered pairs (t.^^,t.) in time, 

b) in the coordinate x, we establish at each partition point 

more than three values of each of the unloiovm 

functions ^lin- ^^.n+l. 

Programs with ordinary and double precision were prepared 
for electronic computers of type YeS and BESM-6 and han be used 
on them without any changes. 
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APPENDIX 


Program \fith Ordinary Accuracy 


Basic Program 


dimension 

lTS = o 
-ITER»7 . 
!TAD=o 
I 5 = 5 o 

T AU= .001 
TKOH=o . 1 

S?S=0.001 

xo=2.s . 

MaS 
N = 21 

DO 2 j=1 


yo<21 f 10> 


OF POOR QDAlIIlO 


u U •• i — * # w 

2 ( 3 .- 54 i 59 *(X 0 * (X1^X0 )‘Fl^AT (I. 

CAUC'HpX^/H,iBjtT«,U,Jt.T E K.. I T,R.., E P_S . T XO T 




. Calculation of “the solution of a 

yoCj;4>5EiP(„9.S696*TKOH>i»Slfii<3 

«f float (1-1 )/ Float <N*'i )) > 

11 yo < I , B> =TKOH*yo ( 1 , A ) 

a . , — ^ - 

S FORMAT. Q 2 H ^comDarison „of _^e 

*'): 0 H Volutions for BP.= TKOH)S 

1 FOR^AT.U) ■ ■ '• 

.. oo:<o'i“i,N 

«0 pRlNj 9 , (yO ( I , J } , Jal , A) 

9 FoRMAy <AE1S.5> 


TAatyfl.,XO^Xl> ^ . . 

control example at moment TKOH 

«j 41 59* ( xo* < Xi -xo ) + 


obtained and.K^ctl 




Subroutine PF 


subroutine FF(MiWN> 
common /HiFy BPjH.TAui, 

1XKC2i),F<iO> , yo<21 ,10> .y<2l ,10^ »yi ^21 *10) 

. NJJ.?, .G.0X0, 1_ 

- Boundary conditions on the right end 

F ^ 1 ' ay V riN , Ji / ; 

F (2Jcy£MNf2) I 




5 

1 f F t N ^ V . G ^ - 

C. Boundary conditions 
> < I / . *y ' 1 1 7 j . 

F {2} a«3 , *y 1 1 , 2 ) +4 , 

5 Right-hand sides of 


on the left end 
*yC 2 » 2 ).,y( 3 i 2 >l 
equations 5 


^ >c( (y <uo1 , 1 5 -2 . *y CN/ 1 > *y CN-1 ' 1) ) y H2+y (N' 2j ) 

**C.55^<(yolN'«'1/l)-2'*yp<NM>-»yO<N-lM)>^H2*yoCN'2>>*0*AS 

3 continue j 
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Subroutine HOK 


subroutine H 0 K{N-Hr is* ITaU' :T 5S'lTa»EPSfT<OH,TAU,yoO'^0"'1^ 
OIMEKsION >' 0 ai 21 , 10 > 

/COHMOn /»!?■/ BP*K(TAUi, 

, 1 XH < 2 i ) * f/<i 0 > * yo <21 ' 1 0 ? <21 Vfl 0 ^ 'yi ' 2 l r 1 05 

T FORHat</> i 

2 f§R^AT<3K n=,I3»3H Ms, 13 , 4(1 1TAU®»-I3' 

, t 6 H ITE'R=* SH- irB«» l3> 

L.. PRiRT, 2 'N»Hf Is, ITaU* ITERi ITB 

3 Format (Sh ■ £Ps- ' El 2' TKOH®* El2,^»5H tali’’' E l 2. A* 
1 AHA 0 =rEl 2 *A,AH)/ 1 “»El 2 .Aj 

prints *£PS.TKOH,hU*XO,Xl 

PR.INT 1 

Ja®® 

BP-0 

/R=(Xi-XO> / <N ^1 > 

/XH<T)=X0 

/ 00 24 I«2,« ■ 

/ IFCTAU.NH.O) TAU1s1./tAU 
IFtTAU.EQ.O> TAU<i»0. 

DO |9 1=1 ,N 

DO 29 J =1 , rt 

23 yo < I , J > =yoo < I , J } 
yi < I , j * 'f > 

29 y<I»J)=yO(** J) 

CALL 0 UT<M*N> . . 

30 IT=-1 , 

Rp = P ■ ' “ ■ * ’ ' , 

I F Up . GT . Ti^ 0 N*i • 00001 ) Goto 72 ■ 

31 IT=IT+1, 

CALC“iSeL<^<M;EPS/A.5 

- r r J Y‘- ^ ^ 

* ' One compulsory iteration 

32 IFCITEQ 1 IERI GOTO'50 
DO |8‘T=i,N 

DO 28 J=1 M 

2S ya*J)=yi(T<J) 

IF(ItB;;^E 0.1 j GOTO SS 

Comparing the convergence of iterations 
■33 00 2 7 ^*1 < 

r F CABS <xH ' J • LE . 1 • E"1 0) GOTO 27 
IFCABsCtytA/jl-y'iClfJJJ/yd^jW'EP®^ 22,27'32 
27 COI^TinUE 

ICCTA-U.e.a 0). fi.eiTo - r- - 

Preparation for the next step 

DO 26 1 = 1 

00 26 J*1 ,1^ 

yo c i , j > =.yi c i , j > 

• ^*'tFtAes<8P-fi<iiH).L£. <TK0H*1 .E-5>5 COTp 58 
25 jB = 0B*>1 

1 F C I T AU . Efl • 0 j GOTO ^3 , , 

IF CIT EQ. 1 ; A.ND,..TA_U„..LT-..t&.QH./1-n^ ) GflTO 41 

It Increase _in .s_tep by TAU • 

^ ^ ■ Decrease in step by TAU 
43 IF<IB.«C. J“'’ GOTO 3« 

Delivery of a result with period IB i 

j B = 0 

CALL 0UT(M*N> 

rtOTO 30 ^ „ 

72 ?AUaTxOH-BP'*TAU 

3P=T<oK 
TAU1 Si ,/TAU 

gOr2 3l 

A1- TAU-1 , 5*TA0 
42 TAU=TAU/2, 
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TAU’l = 1 . /TAM 

If <TAU~TtCO»/iOOOoO.> 

(,ii I . . , . 

■ 6 f 6 r^IAt\( 34H stej? TAU is, smaller -than tiha^peig alffed time 
50 IF (Tau^ EQ, 0> goto 52 

DO 5-i i-i,K 
DO 51 ,1=1 

5’ UWiVo'i'^' 

*9 FOKJ-IAT (39H IT = ITER Iterations do not converse with accuracy _E^) 

- - ., -- - _ - - .- 

PRlJ^T 56, STER 
CAll 

PRXNt 56, sT 

DO 53 i»1 #5 

DO |S 

6 ai.{. out<A/n 5 

STOP , I 

FORMAT (lOH iteration =, ) . 

QAt.1. oVTCMf fi} 

Q 0 T 0 3 1 ,■ 

s? PRINT 56, |T 
56 CONTjfjUE 

CA0& oOr^HtNl 

DO 59 I = 1 , ' 

00 59 J=1 

59 yoo<i,J> 5 »yi <i.j> . 

■ return 

• SNe ■ ■ ...... 


5S 


5 5 ' 
.55 


ORIGINAL I-AGt^ 
OF POOR QUALTT^ 


Subroutine ISCL 


1 02 
1 1 Q 

1 1 1 
161 


1 L 


liHca?) JtiSi ,yoi2il;o] ><21 .io> .vi'2i ,ioi 
HSs'^N , • 

MZo«N"1 

,^f^Zl1cM2*1 . 

K2205S<-MZ 
fil 2 g 1 B 2 * M Z 9 1 
iii?,5ys3*KZ 

MZ31 53’‘M2<-1 

fjO NS «a ^ 

..Ajat.CL A OS 

Cycle in N . , 

jPCri)i1O,lOA,105 
DO 111 laMili/, MZjd 
DO "'ll J = 

62 ( I , J > =0 / 

d|^112^L-1 ■/'2 

Search for. the Pivot Element 1 

I?<A^<|jn'^U?.tT'^PS ccsni)) <i0 TQ 115 
K 5 1 " / 

t» ^W.7 aitllfi ^ 

Pivot element is close to , zero 

■ V E a.b'tCJfii1 f S « 'A *0.^ i JtOtA. • „ 

Rearrangement of the rows and elxmrnatxon 

:*o T'iA^,l = uMZ3i 


:. 2 < 

■ ‘i A 


^Q,L> 50 TS nifc 



116 

1 U 


DO lltj^Ial fM.720 

Go TO 115 
J 2 < 7 >«BZ<bL} ' 

L 1 “'- + 1 

DO 1l 6 Jgll iMZ31 • 

a ,i < { , j ) ^ 3 2 a , j ) ^ B z < 1 . , j J * fl £ 4 ) 

COM I 

A M ? >1 n_-_ 

Formation of the matrices P, 

'ccr :7 To Is: '1 trrs 

XT 5'>2 » < I r Mzsi ? 

JpCIj.cQ,^.; J (JOYO 110 


Q, R 


w* ,J;;37;CI , JV«Z) 

iiQ.o) -eoro ii 7 ' 


117 

11 « 


121 

120 


C 

C 


106 

122 

123 

105 

5*1 


P 2 4.^-1 ..i 

I F 4 M , Li 0 , ^ . . . 

0Z4,i: . «!‘'Z4Z < «J+HS2o)' 

COi. . .-.’eic ■ / 

CC)M 7 iMU£ ‘ 

IF 4 0 iiQ, ~-i 50Y0 106 , 

DO 1 2C t- *=1 ' ('2 ' 

S244,KZ:i1)«Dza*M2f«a51) 

op 121 -;al /hv ,;0 

F, 24i. 

^^0 *2o J~(4£2l#iiZ30 / 

024i,j>»0 

N ».':"1 • ■ , 

End of the cycle in N ■, i ' 

,s 4 n n ^ 0 =, *• -- 

Receipt of the solution * 

00 TZ2 JsI ,ftz - 

00 122 Z=1'M2 

00 123 K»2<‘.2 
00 1 23 .) al /Kz 
D0 l23l = ’i'M2 

'J>3^1 4)4 + 1 , J>'^P24l4,j , i>*yi <K» 1 ) -02 <!4“.1 .J *l)*yi 4 k-1 

Block for calculating coefficients at interior noint of 

calu (:Ftw + i ,nN) / ■ - 1 J. u oj. 

00 5l K'=1,fl2 / / , 

F 1 < ’>' 5 B F ( IL 5 ■ / 

D 0 5 j. I « 1 f. 2 

■?z^bmv/nhin\>dW 4 '!ihnih 

iF(t>y.GG.o> oyacps ■ 

^i^^t2iP=y5N + 2, I> + 0/ 

c A 2 U f f ( K * 1 , JL « > 

y{N?2, X) ay4!’*2, 

0 0.?3 K -f 1 2 1 V» K 4 2 0 


. 1 > 

N 


S3 B2 '■F., I >nCp1 - F 4|>K2) } /py 

n'/»4,'.BS4y <N.>i , ) <.,\ag<y 4iiv2 , • > r 

♦ ACS 4y <N , I ) J +/,;,3 <yD4r.'+i I I J ) J 
IfU>y.EQ,0> 


? + 


F* t** ** 


y 4 N+i , Dv-oy 



ly- 4 y <N . n > +An 5 4 i'<w+i ; -z 
A 3 S < 9 i M-l , r J J « AD^*'<yo < N , I ) > j 
:F 42 y.>^-:.> oy«£?s 


F 4K..y,7.> > / oy 

* 2 ps 


C A t U F F \ ti V 1 ( lU'i 5 f - 

y<N^ lj'=y :>-Dy 

BZ 4 K,HZ 3 l> 6 yi <K-H 2 > ‘ 0 

00 5s IB1 ,117. ^ 

B24A,mZ5i >>2 <K »H231 ) *02 <Kf J)*-y<N + 2# 1> •' •. 
I + bZU, I + M2f *y<N+i ,I)*B 2 4K / I*M 220 )*y <N,I> 

DO 56 K = M7.1 1 ,MZ20 
B255.:'A"3?.*':4,i<>+rAU1 

, Tm. “ 1 >=" T iJC. V tl.^yci t«*4 !-''"!’S2f ^TaUI '• 

Normalization of the rovr ' 




I' f-y. 


ti ^ 1 1 M4 ■ 


0 U . 


: * VI * * T 


V‘ 


> ..uT,^ss<cg(iy ;7 Sot® 
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5 ® 


57 


COfJTlijUE 

If «C 2 c 1 J .^3.0.5 goto IOG 
DO . 

BZ<I, • ■••.'_ 

Mock for calGulating coefficients at the right end 

103 CAtL .-fTTTTN ^KN7 

••■ '■■ ■ DO f (^*>1 ,>^2 

31 " ?'•( <k>' '■■ 

■ DO 3?' I»;i ,nz 

DyctADS<>'tNN,I>)+ASS<y{HK'-1,I>^* 

♦ ABS^y tfiN-a ' I) >->ABS<yo<NN»I>>>*cPS 
lF<0y,CQ.D? Dy^SPS 
V <NN , j } By { NN ; 1 > *Dy 
, tALL- jrftNfV/MNV 

DO .' 33 

SS DI<fTl 5 «ifT(K^-F-CK^>/Dy 

y 1-) >*abs(/o2?5'N*tv 
■ . ir<Py- £C,o^. Dy'=£pG 

y:<NN-H’<i-}sy (Ntf-1 ,io*oy' 


5 ^ 




DO 3 i‘ KbI , ^ 2 ' . . , , 

* B? < -V I ■( '■■■ >' Z-'J ? , 4 f T s k'J rf i i V 07' 
r>y t= ^ AB s < y C H Nr"2v -I:) '> + A B 3'i Y'i fj'N '•'H Z 
+ A B s ( y ( N fg - 3 ■' 1 0 ■) AB5 ^ y ■< M‘N’"' 2 .y '1 3 ') ■-> 
iF'toy.Ea.o'^' Dy»^ip-s- 

y:</jf <--2 1 1 )=y’<NNK 7 ;,,i:) +d/- 
CA^:U fF<NN’«N'' 


‘ . y^<l<\-’2 * I .) N N'‘-*2^ -I^-Dy- 

D.O. 3 'S' K?--i . . 

C 3 -^'S/ “? <A‘) 

^ 10 ¥h:l-iM.‘:t TVM 2 -Eg')-*/^ 

, . Normalization of the row 

n *0 3 y IBT-J,nz' 

D0^3s’‘StrY,AizVo'- , 
jy:tABs-{.0'=2-C.V»>>=>h'lvTlA''?S‘C(^z:tV>^^ e S5: 

, , c 2.C 1.) :•; j-)' 

3S s^dfo- vo^o-'- 

57 .____•. 1 ... .. 

GoT(y 'BXock for calculatimg coefficients at the left end 
1 c A'lrU* 

, o'o- <i*V K'’-V,h'Z' 

4-1 Vt 

DO *'A' 

D/«^;A'kS (5 , 1) > <.-AffS <y ^4, J-yi 

• +‘AB5-^y^Z . I-).J+'ABS (yft<3< l> ) > 

IF(.Dy.EO.o? .oysEPS 

• y c 3 -f cs ' 1 ) '('Dy 

CA-ki? ,frp(i ,KI^) . 

yc3',':-j--y.<5 » 1 5 fay ,. i 

• , d'O A'3 ,K«M2] 1 » I^ZgO .. _ 

. K| ii ' 

' ly < 1 . i ) >.'fABS <yo <2 > tyy> '^BPS 
iFcay..-EQ.o' .oy^eps, 
y(2»iy^y<?'i-)^Dy 

y.c 2 f ij^y^E' J)'*ay ■ 

00 «2 R«f^ 2 n,M 2 SCi 

tfFgz CK , RZ ) » ( I* 1 .; <» t « Rf w j> $ /fry 

^ADSCy (3 <n>^ADS(y6<1('i>'y)*Bi*i 
lF<Oy,-r:Q, 0 > .&y«2PS 

CABU cF( 1 ,N«> 

■ynfn“ytii i)-fry ..- • 



Lk f>Z(K,t+HZ2Qlatf1CK-«Kj>-»FU-HZ>^/D> 

45 fir 

1 Bi C K . 1 ■* H 2 y ,< 9 >.li*-8 2(,KfI*WZgo5*'^l'‘' 

/Elimination o|LA(0,5_i . 

nO- ^0 J “HZ'* 1 I MZZO 

IF<?Z<J » K> > £q, 0> - epio 

DO I si . 


• 6 . 0 . 


61 

60 


CONJinUE 

DO 6s 1 = 1 , f^Z 

82 ta^j) aBzU 

8Z<j'.I-»MZ5=e2<j, I+MZ20-) 


6Z BZ<4 

^0 ^ r^r^iriT'r-rTTrrzT -■ 

FffAlsJaUi^j?) .uT.abs<cz(U>) 

bO 0.) COTO ,100 

^ , Vatlvite in 'deliveiT,] •‘■' 

foT')>i?iNr 


goto 4S 


DO 3 ‘ 'Wzib 


P8-1.NJ-:-?. 

1 ; FOBMA.T 

2lo>'l^'AfitO£V2?4) 


, ISfLii--) ' r— ni 

i723E Matrix is singular !_P' 


N =,„13aJ 


Subroutine OUT * 


SU8AouTlN£„0UT<M,fj) 


1 foU'iK'fl ' 

5 J Q R M.4Jf-.f .J-l-B-''-'' 5 

I' FOBMAT (lOH: 


Tirae ,il2,4) j 


I 


p-K i'l' X y / p'p 

PRINT 1_ 

io print L U'l > / (^1 < J 
print ^ 

RETU«H 
E«0 ■ 



Program with Double Precision 



0OU3i,e Pr = Cisjo^N / 0 , E P S , T 0 H , T A U , X 0 , XI 
dimension XocE, 1 , lOl 
• ITO =0 / 

ttER=7 /' 

i " A U = 0 j 

I B 5 0 / 

TAO=.oOI / 

t?s= 0 , 00 l / ' 

X0=0 . s • / 

/ • ■ 

M = 2 / 

N = 21 / 

00 ,2 1 = 1 /N, 

^ yo‘*'25=6SlNC3.14l59»CXo + < X -j « XO > * P UOA-T <1 ~-i>/|:L0AT(N-OI) 

2 yo«i,i)=o. ■ - . 

C 4 L Lxi.M.O X M . T.5,^ <■ T 4J.I * T T 5 O . 5 T ^ -S_kJ /• , T C A LL . T JI.I I-. V , . i ' 

Calculation of exact solution of control example at moment TKOH 

J^O^I«‘{,'^=6^?^?<-9,e696*TKOH>*OSI'<<3.lV'i59*<XO*CXl-XOl* 

DRXJviT t . _ . 

8^ FORMAT (32H Comparison of obtained and exact ' 

2 OH solutions at BP = TKOH) 

' r V • I ' / V / - . ■ 

D'b 1 0 I =1 /'N ' ■ ■ ■' 

VO-pPtNT 9. (/OCi, j) , J = i ,4) 

9/pORMaT <4'tl5.5> 

./ EHO ■ / 


Subroutine FF 




-H,TAui »XH,F*yO'X,yi 
common /HIE/ 8P,h,.TAui, 

1XK^2l),F(iO)^yo<21,‘'o>,y<21,ioJ.yi^Hl,10) 

right end > 

F;i>=y<NN'i) . 

goto 3 • ' “ ■ 

1 T C ( N jj£ . «s ^ n T A. -s _ 

Boundary conditions on the left end 

P ^2 j =-3 . *y ll , 2) *.4 , *y <2 , 2 j ^y,{| ^ 2) 

Right-hand side of the equation 

Z COf'iINOS 

FtI)B<ty(N+1,l>-2.*ycN)l)+y(N,1,i)>/H2+y(N(2>) 

5**g,-.^;'-j^^O<^*1.2I-2.*yo<N,2)+yo<N-1,g>>/N2*0.45 

RSTUrn “ 

£N0 ■ 
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Subroutine DHOK. 


'SUf'ROuTiNE DH0 K<m,M / IB / ITftU, iTERf lTB,EPS,TKOH,'rftUi>'00^^O*!^1? 
OOUBUF. P;i£ClsI0N /O 0 , £ PS , TKOh , TAu f , X 1 
double B'RECIsION BP/H,TAu1 ;Xh, Fiyo'^/yi 
OIHEnsIO.X yOo(21 ,105 

common AkiF/ BP,H,TAui, . ' 

JXH<ei)/F(1'>5,y0(2V/l0>,y(21,‘l05»yH2l,10> ' 

i foRmatA/) I 

P R 1 M T / 1 

n=/I3f 3H M=,I3,a IB=/I3<6H 1TAU=,13/ 

16H It/eR=> , IT3=/I35 

PR IM.T2 F N , H' I e , n AH * ITER F S TB ... -- ' 

P R J Nt 1 ^ 

5 FORflAT<SH £Ps=/Di3 .^F , 6 H T X-OK = , ^5 1.3 . f'S H TAU=/0l3.4/i 

1 tHyXo= ' &1 3 . A , 4H Xl =-F D 1 3 , A) v . •'• 

PR / ^13 F EPS”/ T KOH , TAU F X 0 F X 1 ‘ ‘ 

prAnt 1 , ; 

j^=p ' •. ■ ; 

B ^=0 

/XHU)cXO; ' 

/ 00 2 a 

I XH”)=XH<l"1)^K I 

1 F <TaU'/NE . 05 TAU1=1./TAU ' 

IF<fAU/'EQ.u)TAUl= 0 . 

DO 29/1 = 1 
DO' ?9/ J =1 > M 

* y 0 <VJ) = y 0 OU » J> 

X 1 < y / J > = y 0 ^ 1 , J 5 
> y ( I/' j ) ^yo ( I f j ) 

CAl/u dOiJTCM,n) 

) lT /=-1 ^ 

B^psBp<.TAu „ 

JFtSp.GT.TKOH+O. 00001) GOTO 12 
/cALt jllCLXilxA^JF/SjlA . 

,.1 One compulsory iteration t' ‘ . 

: I F ( I V-c-^-i-r crKv~‘'a*Orti '5"0" ' 

DO |a 1= 1 , M 

DO 2 a J = 1 , M 
y <l F j >=yi (I , j ). 

I F (ilB /EQ . 1 > GuTo' 55 - 

„ ! Comparing SJtie convergence of Iterations : 

DO 2 y — “ — ~ — — 

DO 2? J = 1 ,r;i - - ' ' 

lF(OABS<y<I,J) >...lE;T;d~2o5 goto . 37 ... , 

-If|5PABS<(.y<r , J 5-yi n,J))/y<I,j)5-EPS)37/27F3 2 

] F ( 'f A II ■ F Q Cl3 — c.ax.a — ILJ 1 

I Preparation f^r tiie next step 1 

DO 2V f -1 , IN — • ■ — 

DO 26 J =1 ,M 
y 0 C i , J 5 ='/1 Cl , J ) 
y<i j> = >;i FJ) , 

I F CUA 8 & CbP^TkOH) . LE * 0 • 00001 *T KOH ) GOyo 58 
0070 ^3. 

I F < T t_E.b <t_-.Ajo.r 5 .Tjt..u L_t_ T_!CftiF_Aj<^« 1 ,?.n.T.Q. a < 

■; In crease 'in the step by TAU , ' 

lP.<i — : — • — — - — -- 

I Decrease- in the step by TAU _ 1 

1 F ( I B . Nl.._J_D;/;GOLOl_3jDr_._/ 1/ J ^ 

; Delivery of the result with period IB i 

CALI. oOHT<M,N) 

goto 30 ^ • •• ■ 

TA'UFiTxOH-BPtTAU— ’ 

.be-^ikoh-- 

'TAU1=/, . /tao 

goto 31 . 

TAlj = 1 .S*TAU 


29 



TAUJs, . /taw 
ISOTi aS 
.2 TALI = TaU/2, 

TAUIsi./tAu 

I F < r AU~T *0?/ 1 OOOoO > 5 AAfA3f43 

/ T N ^ 5 ... 

6 FORMAT (34H Step TAU is smaller than the permitted time =, , 
n 7' 6/ BP*; T Au" 
s. TAU.EQ.®> goto 52 

DO 5 i 1 = 1 ,N ■ ' 

OO 5i J =1 M 

Evy <i ' j)=yoi?»J> 
goto a2 

52 S>».I NT ,o 

9 FORMAT (39H IT =: ITER Iterations do not converge V7lth precision JIPS) 

PRlNf'‘‘56) iTER ' • 

CaLU DOyTCM.N) 

PRINT f6/lT 

DO 53 1 = 1 , N 

• 00 53 J=1 « 

53 yi <I| J> 

• CaLG 00uT<HrN> 

STOP • • ■ 

55 P.R1.KT 5-fe.TT - - 

56 FORMAT (lOH iteration =, 13) , 

goto 31 ^ } 

5? print 56. iT , I 

58 CONTxmUE 

CALL dOUT{M,N) 

DO 59 I=T.N 
DO 5 9 J = 1 ( 

59 yoo ^ I , J> «y 1 1 1 . j 3 

RETUftn 

END 


Subroutine DSCL 


• C 


SuSKouTINE^OSCUHM, NN, EPS) ■ , ^ „ 

.double precision BP r H , TAUl , XK , F .yo ' . y 1 

COHrtON 7HlF/ BP.H.TAQi, » , . 

TXHt2l)/PnO),y0<21.l0>.yC21.1D>ryl'2'l.‘S0> 

mz = mh 

M2TT=MZ*1 

KZ20_2*HZ 

HZ2T=?*M2<-i 

mZ20=3«M2 ■ 

MZ31s3*M2*'> 

N=fU-i 

a . * rtj! . 

Cycle in N 

162 lF^«)np,io'“,iq5 

110 00 111 IsMAij.MZjQ- 

DO 111 J=1 »Hi3l 

111 BZ < £ f J > =0 ^ 

101 DO ">12 Ld.MZ 

...an.i.c.0 ... 

Search for a pivot element 

I F <dA|s ( bM'i ‘ t > > . Lt • D^BS (C z ( 1 ) > > 11^ 

K = I 

C2<1>s.BZn,l> 

^ /\jj T ♦ 4i IJ P » , 

Pivot element is close to zero 

IF^UABSICS'Di.VE.I •O'-SO) QWTU 10 0 


1 1 
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.■Rearrangement 'of the rows and eliiiiihatlnn 
— D'cr- T i -*r-rsTrpu'3T 

a “T f tC t \ , f 1 \ . 


“D'cr- T i *r-rsTrpu'3T 

BZCK, J> = 82U,J) /C2<1) 
IFCi^lEQ, L> ao T0-'i1<> 
-CZC2)«BZ{KJ;)> 

BZCK,j> = bz 51 .,J) 

B2U, j> = CZ<2) 


Dfc ^ ^ t yj ' Vi 

C0nT2mU£ 
BO 1 •! 5 l~’ 


C2 Cl >'«=92 a f L> 


ns 


11^ 
in 
1 12, 


118 


121 

120 


DO 1lA J = LnMz31 I 

32U , J>=BZCI ,5>.BZ<l./g>*C2Cl> 
continue 

c,oNTtm.US , 

Formation of the matrices P, Q, R 

-D-o— rro-Tr=rrT<Tz- ^ 

IF CN. EQ. -n GOTO 118 
DO M? Js1/H2 

P2 CN-s-i f I , J nsZ ( I , J'>Mz > j 

IF(N.E«.0)„6otO 11T ■ • 

QZ<N, J n>B8Z(j,g4.MZ2o) 

continue 

continue 

IF<N,EQ.~1 > GOTO 106 
OQJ20 T = nBZ 
BZ(r,ji,Z3l)=8za<-MZ^M23l>‘ 

00 121 J«1 'H220 

BZn,j>=BZCl->MZ, j*H2> 

00 1 20 J=M2£i ,M 230 

8Z<I,J}=0 

NsN'^1 %■ 

Sid of the cyc le in N 


ORIGINAL PAGE Ig 
OF POOR QUALIT5? 


F » 


.Receipt of, the solution 
D0 'T22 JV1 ; M2 
00l22nt=1'Mz 

DO 123 J=1»HZ 
DO 123 2 = n«2 

'4’^ V -J ZJ L- ' . J .t) • ^ H K - 1 , ! ) 

at >te^£or polri-'ol M 

■00 Sl K=1 ,M2 
5l -FI <K)=F< k> ■ 

OO ?2 1=1 ,MZ 

oy“CoABS(ylN + 2,n)+DABS{/(N<.5,I)5+ ' ' 

'*iF coyl^^Eoto? 

ycNn p=y£^^2,I>-0>' 

00 S 3 K=M2l1 ,M52 o ■ ' ■ 

<hCK-MZ)-l 


1 06 
122 

1 23 
105 


53 BZCiC 


54 


52 


5S 


lF<Oy.QCl.o> Dy=EpS**2 

l?« i"' 

t) + oy 

ycNf n=y?ni i>~oy 

DO ?2 K=MZl1,MZ2o 


I) 


31 . 



00 So f = M2H,MS,20 

56 B 2 < ‘= 5-7 c,v:>.mj 731 ^ i.yo.f ii-M *T/\yl 

I', Normalization of the row , 

DO S/y — ^ff^TnrTTTiTSYD .r.’---- - 

C-ZWOsO ■ • / 

DO/SS J<^1|t^230 - - - /^» . . 

IF'^OAbS(32VI,J)5 .LT, dabs; CZ (1)3 > 59^0 '58 
C/ 2<1)=8Z5 j/J) 

SB.CONTjNyB / 

/if (tZ < t > /eQ. 0 , > flOTO 100 
' / 00 ?7 J=/i ,MZ31 

^/7. BZ(i ,j^=Bz<.i, j)/czn> 

; B'lockj;^_ 5 al< 5 ilat^^coem 

io'3 cAL'-/Tr(W'irN) ' - -•''••:•> . 

DO 3-, K=1 ,M2 

gl Fl(i'>aF(K) 

.DO ^2 . 

^ 0)'*‘0A8S(>-(KM,X))+0A8S<7fNM-i,t?)*, 
<*DAB5.(>'(NN>-2»2>>'*’OABS<yo<NN»n.)3*EP5 
irSCPy.Eo.o) ovsEps^^a ■• ■■ 

y (MN, I ) =y <hN, n ^-dX 

CAUU FF(NNfNN) 
y CN«, i>=y<M^, nr-D/ 

00 33 tC = l ,«2 , . — 

55 8Z(^,]) = <F'i (<)~F(K) 3 /OX- ' ■ 

Dy=(0AS5(X''^N-1.I))'*-DABS<y(UN,i))4 
■ *DAB5(y<NN-2iI)>*DABS<yo(NN-1,l33)«HPs 
jpCOy EQ.o 3 oysEPS^^S 
y (NK^i , n = y CNN-1 , D + oy 

CAl-l- ffCNNiMN) 
y C.MK-, ( DaX <NN-1 , I) "DV 
DO 3a Kai,MZ 

34 B2 (F- , S-i’WZ) = C FI (IC) -F < K) ) /DV 

0 /= dabs cy CNN-2 ,I>)'«-DABS(yCNN-T'l)3 + 

♦ DABC’ Cy <NN-3* .))i-OABS<yoCNN'-2fI>3 5*fcPs 
IF(Oy.EQ.03 Dy=EpS“*2 
y (NN-2 I n =y <NN-2 , 1 ) +DX , •■ - 

CAU. fF<NNjNN> -• ' 

y (NN-g , I ) =X CN^^-2 , I ) ”0y ■ 
nO 52 / "'Z 

32 82C*^,I+M220) = (F1 (iC)-F<iC)5/Dy 
DO J S Kd , M2 
BZCA MZ31 >=F1 CK) 

DO^gla', ,«Z , •• 

55 B2^A,nZ31 >=B2(IC.M23, >*82CK,I><'X<NN(U • ■ ; 

1 + 8 2 C ^:jX*ji.Zi.s2CC.!J 4 / 5 - 1 . 1.5-tB.Z.rK-^ 5 X ( N N • 2 » J ) 

Normalization of the row , 

no 5 7~l»r,'i'r« 

CZ<1)‘S0 

DO 5b JbIjMZso 

IF COASS(32Cj , j)} , lT,DABS<C2Ci>)3 c 9'^'0 33 
cznjcBzcifj) 

38 CONTjnUE 

IF <<^2 (I > .5Q.0. ) flOTO 100 

00 5? 0 = 1 tf^231 

57 eZ<d J>:S2^I . <>>/'CZ(l) 

i_Block for calmlating coeffic^nts^ on the left end 
10A pAU^ ' - , .- 

00,51 

t1 . FI (K>^FCK) 

Dn4Ai=1>«Z' 

Dy=^DABscyc3, ni + oABscycg,!))-!- 
«.DA0S <y C A , I > ) +dABs O'O ( 3 f 1 ) ) > *EPS 

1 F CPy EQ . 0) Dy=EPS*‘'2 • 

y<3 ,2)=y<3, i)*oy 

CAtU pF<1 ,NW> 

DO fc3 K=M21 1 ;MZ20 

43 B2<>i/l) = CF1(K-H2>-FCK-M4>)/py 


■I 


ByS‘'BA3®<Xj?'’^^i*0AB3Cy<3iJ))^-- 
*DABeC?Cl,i)) -v-dABs <70( 2 i I > > > *eP® 



lF<Oy.5Q.o) DyaEPS*»2 
y <Hi n=y <2Fn»D’/ 

CAli. ?P<wWK) 



/ (2 1 n'ay (a f I > -oy 


A2 



. cy(j(< 

I F '■py , 20 . 0? oya£ps*«2 
y(j> ' s>voy . .. - . ■ - 

CAl Fr < -i f NN > . . 

yM ; n>=ydi d )-[>/ - . . 

r »0 '‘A K=>'iZ 1 1 • ■ ■- ■ 

50'»SKart2idHZ2o 

>’=?1 <K"M2) 

DO A 5 I = 1 . nz 

Ellmin^ion of A(0) 

DO 4>0 ■•j’ = M2’i 1 ,HZ2o 

rFd;a,<g,K5 ;eq.o} goto 6o- ■ 

00 «■} ia 1 ,f'i 2 
BZd', ^ + 

BZ^'J ,MZ3i > =S2 ( j ,H231 ) <3 »K1 iiB2< J ,K> 

CONTjjjUB 
DO «2 a«1,N2 

&0 » MZ2o 

sz , : >«3k u , 

SZ-JJ , !-4!-l7.):62 < J ( I+M220) 

g 2 W f ^ _ _ - -■ 

Normaliztion of the row - 

pO X"fT£iT,K2'20 

b 2 < 1 > = 0 . 

DO'*3J=‘1,'^220 , 

iFdA3s<^2drjn.LT-o^ss(c;<u>> g^^To as 
J) 

CO'^TlNliS 
rFCCzid 







:Q . 


0 , > GOTO 100 


1 £>0 

3 

1 

.2 


00 Ay 0 = 1 ,M 23 i 

Failure in Delivery ;- 

PRINT • ;« , w . K ■ 

-00 3 ;=1 ,m2So 

— i^./. /-- . -J.A . — - 

FORMAT C/23E Matrix is s ingulf: W 

- an “•« , i J / a "' R -^ jAS ^ H-n aAT =) 

F0!^MAT<lQi)l2.A> 

STCP 
HNO "■ 


= ,13. 


Subroutine DOUT 


subroutine 00UTC«,N? , ^ . V, w 

double precision BPjKfTAUI f XH. F'yo'y.yi 

- COMMON /HI?/ BP,H,TAui, 

1XH<2l> .FCl-OJ ,yo-<-21 d0> iy<21 ,10> 'yi^Zl dO) 

1 FJORMAT^y/J ■ 

^■2" F0«flATdlDd .3> ■ ^ 

9 FORf^ATdoH BP£Hf5=f 01 5, 
print 9, op 
PRINT 1 
00 1 0 I»1 / N 

10 print 2,X’rld> d/l «I' J5 Wod«> 

1 , 

Rc'i'URW 

ENO" 




Printed at the ICR of the AS of the USSR' 
Sent to press 3/20/78 i 
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